NUMERICAL COMPUTATIONS WITH i?(div)-FINITE ELEMENTS 
FOR THE BRINKMAN PROBLEM 



JUHO KONNO* AND ROLF STENBERGt 

Abstract. The i?(div)-conforming approach for the Brinkman equation is studied numerically, 
verifying the theoretical a priori and a posteriori analysis in |27II28| . Furthermore, the results are ex- 
tended to cover a non-constant permeability. A hybridization technique for the problem is presented, 
complete with a convergence analysis and numerical verification. Finally, the numerical convergence 
studies are complemented with numerical examples of applications to domain decomposition and 
adaptive mesh refinement. 

1. Introduction. The Brinkman equation describes the flow of a viscous fluid 
in a highly porous medium. Mathematically the model is a parameter-dependent 
combination of the Darcy and Stokes models. For a derivation of and details on 
the Brinkman equations we refer to |29| [111213 |3S]. Typical applications of the 
model lie in subsurface flow problems, along with some special applications, such as 
heat pipes and composite manufacturing [26^ .20^ . The effects of taking the viscosity 
into account are most pronounced in the presence of large crack or vugs, typical of 
e.g. real-life oil reservoirs. The Brinkman model is also used as a coupling layer 
between a free surface flow and a porous Darcy flow |17| . Numerical results for the 
Brinkman flow have been previously presented for the Hsieh-Clough- Tocher element 
in [38 , for the Crouzeix-Raviart element in [5^, for the Stokes-based elements with 
various stabilizations, including interior penalty methods, in [2T1 [T51 [T^ [TUl [TT] . 
and for coupling the Stokes and Darcy flows with an SIPG method in |5S]. For the 
iJ(div)-conforming approximation, numerical results with a subgrid algorithm can be 
found in |;23|, and with modified i?(div)-elements in ^371131]. 

The structure of the paper is as follows. In Chapter|2]we briefly recall the mathe- 
matical formulation of the model, and introduce the necessary function spaces. Chap- 
ter [3] carries on to introducing the i/(div)-conforming finite element discretization for 
the problem, along with the Nitsche formulation for assuring conformity and stability 
in the discrete spaces. We also recall the main results of the a priori and a posteriori 
analysis carried out in |27| . along with the postprocessing procedure necessary for 
the optimal convergence results. The results are extended to cover a non-constant 
permeability. In Chapter [4] we introduce a hybridization technique for the parameter 
dependent problem based on previous hybridization techniques for mixed and DC 
methods [TBI [HI [5] . The practicability of the hybridization and the benefits therein 
are discussed briefly. 

We end the paper with extensive numerical tests in Chapter [sj We flrst demon- 
strate the convergence rates predicted by the theory for both the relative error as well 
as the a posteriori indicator. Furthermore, the performance of the method is com- 
pared with that of a MINI flnite element discretization. Next, the importance of the 
postprocessing method is clarified and convergence of the hybridized method is stud- 
ied. We also apply the hybridization procedure to domain decomposition. The weak 
enforcing of the boundary conditions and adaptive refinement techniques are studied 
in the framework of a flow in a channel with a parameter-dependent boundary con- 
dition. The chapter ends with a practical example of the Brinkman flow with actual 
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material parameters from the SPEIO dataset |13j demonstrating the applicability of 
the estimator to adaptive mesh refinement. 

2. The Brinkman model. Let il C M", with n = 2,3, be a domain with a 
polygonal or polyhedral boundary. We denote by u the velocity field of the fluid and 
by p the pore pressure. Let K denote the symmetric permeability tensor and /i and fl 
are the dynamic and effective viscosities of the fluid, respectively. With this notation 
the problem reads |35l \M\ 

-flAu + ^iK^^u + Vp = f, infl, (2.1) 
div u — g, in il. (2-2) 

To simplify the mathematical analysis we assume the permeability to be of the 
following diagonal form 

K-^{x) = a{xfl, xen, (2.3) 

in which cr is a strictly positive, piecewise constant function. We furthermore assume 
both the viscosities to be constant over the whole domain CI. By scaling the equations 
with the dynamic viscosity, we arrive at the following scaled problem 

-t'^Au + a^u + Vp = f, in fi, (2.4) 
div u = g, in il. (2-5) 

Here the parameter t represents the effective viscosity of the fluid, whereas cr reflects 
the variations in the magnitude of the permeability fleld. For simplicity, we con- 
sider homogenous Dirichlet boundary conditions for the velocity fleld. For t > the 
boundary conditions are 

u = 0. (2.6) 

For the limiting case i = we assume the boundary condition 

u-n = 0. (2.7) 

For t > Q, the equations are formally a Stokes problem. The solution {u,p) is 
sought in V X Q = [H^{n)]'^ x LQ{n). For the case t = we get the Darcy problem, 
and accordingly the solution space can be chosen as V x Q = H{div,il) x Lq{CI) or 
V X Q = [L^(ri)]" X [H^{fl) n Ll{Q)]. Here, we focus on the flrst choice of spaces. 

In the following, we denote by (• , • )d the standard L^-inner product over a set 
D C M". D ~ fl, the subscript is dropped for convenience. Similarly, (• , • )b is the 
L^-inner product over an {n — l)-dimensional subset B (Z Cl. We deflne the following 
bilinear forms 

a[u,v)=t^[Vu,Vv) + {a'^u,v), (2.8) 
b{v,p) = ^{dw v,p), (2.9) 

and 

B{u,p; v,q) = a{u, v) + b{v,p) + b{u, q). (2-10) 

The weak formulation of the Brinkman problem then reads: Find {u,p) E V x Q 
such that 

B{u,p;v,q) = if,v)~{g,q), y{v,q)eVxQ. (2.11) 
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3. Solution by mixed finite elements. Let ICh be a shape-regular partition 
of il into simplices. As usual, the diameter of an element K is denoted by hj^, and 
the global mesh size h is defined as /i = maxKeKh ^k- We denote by £h the set of all 
faces of ICfi- We write He for the diameter of a face E. 

We introduce the jump and average of a piecewise smooth scalar function / as 
follows. Let E = dK D dK' be an interior face shared by two elements K and K' . 
Then the jump of / over E is defined by 

m^fw-fW'. (3.1) 

and the average as 

m = \{f\K + f\K'). (3.2) 

For vector valued functions, we define the jumps and averages analogously. Denoting 
by n the normal vector of a face E, we define the tangential component on each face 
E as 

Ut = u — [u- n)n. (3.3) 

In addition, we assume that the piecewise constant permeability field agrees with 
the triangulation JCh and that there exist a constant C > 1 such that 

^<^<C\ Vif,i^'e/C,. (3.4) 

On each edge -E € £/i we define — (cP'\k + o'^|a'')/2- Thus we have cP' ^ cP'\k and 
(T^ ^ cr^lx' for an arbitrary face E. 

3.1. The mixed method and the norms. Mixed finite element discretization 
of the problem is based on finite element spaces Vh x Qh C i?(div, O) x L\(yi) of 
piecewise polynomial functions with respect to ¥ih- We will focus here on the Raviart- 
Thomas (RT) and Brezzi-Douglas-Marini (BDM) families of elements [8]. In three 
dimensions the counterparts are the Nedelec elements [32] and the BDDF elements • 
That is, for an approximation of order fc > 1, the fiux space is taken as one of the 
following two spaces 

V^^ = {v^ H{div,n) I v\k e [Pk-i{K)r®xPk-i{K) VK e /C,,}, (3.5) 
V^^^*^ = {ve H{diy,n) I v\k e [Pk{K)r yK e ICh}, (3.6) 

where Pk^i{K) denotes the homogeneous polynomials of degree k — 1. The pressure 
is approximated in the same space for both choices of the velocity space, namely 

Qh = {q^ I q\K e Pk-i{K) -iK e JCh}. (3.7) 

Notice that C Vj^^'^^ . The combination of spaces satisfies the following equi- 

librium property: 

div Vh C Qh. (3.8) 

To assure the conformity and stability of the approximation, we use the an SIPG 
method |24t[33|. also known as Nitsche's method, with a suitably chosen stabilization 
parameter a > 0. We define the following mesh-dependent bilinear form 



Bh{u,p; V, q) = ah{u, v) + b{v,p) + b{u, q), 
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(3.9) 



in which 



(3.10) 



^ {^(Ki, KDs - ({|^|}, KDb - KDb} 



9n dn 

Then the discrete problem is to find G Vh and ph G Q/i such that 
Bh{uh,Ph]V,q) = (/,t>) - {g,q), V{v,q) e Vh x Qh- 



(3.11) 



We introduce the following mesh-dependent norms for the problem. For the ve- 
locity we use 



KeKh 



E iiv^k+Et^ 



EeSh 



and for the pressure 



\\a-,t,h — 



/)2 



(72 /i2 + 



i^pwIk + E 



2 

O.E- 



(3.12) 



(3.13) 



Note that both of the norms are also parameter dependent. To show continuity, we 
use the somewhat stronger norm 



EeSh 



du 
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(3.14) 



It is easily shown that the norms (3.12) and (3.141 are equivalent in Vh- We have the 
result laSj, with C/ > . 



hE\\g^\\lE<Ci\\Vv\\lj,, yveVh 



(3.15) 



3.2. A priori analysis. For the proofs of the following results, see [27]. First 
we note that the method is consistent. 

Theorem 3.1. The exact solution {u,p) £V y.Q satisfies 



Bh{u,p; V, q) = (/, v) - [g, g), V(t>, g) G F/i x Qh- 



(3.16) 



In addition, the bilinear form a;i(-,-) is coercive in Vh in the mesh-dependent 



norm (3.12 1 



Lemma 3.2. Let Ci he the constant from the inequality (3.151. For a > Cj/A 
there exists a positive constant C such that 



ahiv,v)>C\\v\\l,^^, yveVh. 



(3.17) 



With a trivial modification of the proof presented in [27], we have the discrete 
Brezzi-Babuska stability condition. 
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Lemma 3.3. There exists a positive constant C independent of the parameters t 
and a such that 

sup > Cl\q\l,t,h, V<? e Qh. (3.18) 

By the above stability results for a/i(- , • ) and b{- ,■), the following full stability 
result holds, see e.g. |B]. 

Lemma 3.4. There is a positive constant C such that 

sup II f''^'''f'r'if > C{\\r\\,.t,k + \ls\U.t,h), V(r, s) e V,, x Q^,. (3.19) 
{i!,q)eVhXQh \\'"\\'y,t,h + \\m\a,t,h 

In iJ(div), a special interpolation operator Rh : -ff (div, il) Vt is 

required, see |S]. We denote by Ph : L'^(57) — > Qh the L^-projection. The interpolants 
possess the following properties: 

(div {v - Rhv), q) =0, Vg e Qft, (3.20) 

(divt;,<z-P^q) = 0, Vi^en, (3.21) 

and 

div Rh = Phdiv. (3.22) 

By stability and consistency we have the following quasioptimal a priori result shown 
in [27]. 

Theorem 3.5. There is a positive constant C such that 

\\U - Uk\\a,t,h + jPhP - Phla.tji < C\\U - RhU\\„,tM. (3.23) 

This contains a superconvergence result for — Php\lcr,t.h, which implies that the 
pressure solution can be improved by local postprocessing. Given full regularity, we 
conclude the section with the following a priori estimate. 

Theorem 3.6. Assuming u e [i/'=+i(r2)]" or m e [i7'=(f])]" for BDM and RT 
elements of order k, respectively, we have 



C{ah + t)h''-^\u\\k, for RT, 

C{ah + t)h^\\u\\k+i, for BDM. 



3.3. Postprocessing method. We recall the postprocessing method proposed 
in 127] based on the ideas of [SO]. Due to the varying permeability parameter a, we 
modify the method accordingly. We seek the postprocessed pressure in an augmented 
space Q'^ D Qh, defined as 



, ^,{q€Ll{n)\q\KePk{K)yKelCh}, forRT, 
^{qeLl{n)\q\KePk+i{K)yKeICh}, for BDM. 



The postprosessing method is: Find G such that 

PhPl^Ph (3.26) 

{VplVq)K^{t^Auh-a^Uh + f,Vq)K, ^q ^ {I - Ph)Ql\K- (3.27) 
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In [27] the analysis of the postprocessing method is performed by treating it as an 
integral part of the problem by embedding it into the bilinear form. Note, that this 
is solely for mathematical purposes, in computations the postprocessing is performed 
after the solution of the original system elementwise. The modified bilinear form now 
reads 

KeKh ^ ^ 

(3.28) 

The postprocessed problem is then: Find {uh,pD £ Vh x Ql^ such that for every pair 
{v,q*) e Vh X it holds 

BUuh,Pl;v,q*) = Png; v,q*), (3.29) 

in which 

CHif,g;v,q*)^(f,v)-ig,q*)+ ^ ^/f"" {f ,V (I - Ph)q*)K ■ (3.30) 

Now using exactly the same arguments as in |27] for the case cr = 1, we can show 
that the solutions of the postprocessing procedure and the modified problem (3.291 
agree, and we have the following quasioptimal a priori result. 

Theorem 3.7. For the postprocessed solution {uh,p'^) it holds 

\\u - u^M^,t,h + l\p - P*hL,t.h < C inf {\\u-Rhu\\,^t.h + \lp~q*l\a.t.h (3.31) 



Assuming full regularity, we have the following optimal a priori result for the post- 
processed problem. 

Theorem 3.8. Let us assume {u,p) e [H''+^{n)]" x H''+^{n) oruG [iJ'=(17)]"x 
xij'^+^(51) for EDM and RT elements of order k, respectively. Then we have for the 
solution (ufnP'^) of the postprocessed problem (3.291 

uk+l 

\\u ~ Uh\\a.t,h + \Ip - P*hL,t.h < C{{<jh + t)h''-^\\u\\k + -^—\\p\\k+i}, for RT, 

all + t 

and 

\\u - Uh\\a.t,h + \Ip - P*h\la,t,h < C{{ah + t)h''\\u\\k+i + ——\\p\\k+2}, for BDM. 

an + t 
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3.4. A posteriori estimates. In this section we introduce a residual-based 
a posteriori estimator for the postprocessed solution. It should be noted that the 
postprocessing is vital for a properly functioning estimator. We divide the estimator 
into two distinct parts, one defined over the elements and one over the edges of the 
mesh. The elementwise and edgewise estimators are defined as 

/j2 

Vl = ^2j;2^\\ - t^^Uh + a^Uh + Vpl - fWl^i, 

+ {e + a^hl)\\g-Pn9\\lj, (3.32) 
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|2 

\0,E 



MM 



0,E- 



The global estimator is 



V = 



HE ^?.+ E 



1/2 



(3.33) 



(3.34) 



\KeiCh 



Ee£h 



Note that setting t ~ gives the estimator of PO] for the Darcy problem. In the 
following, we address the reliability and efficiency of the estimator and show the terms 
of the estimator to be properly matched to one another. The estimator introduced 
is both an upper and a lower bound for the actual error as shown by the following 
results, provided that a saturation assumption holds. The arguments presented in p7| 
hold step-by-step for the estimators (3.32) and (3.33) in conjunction with the scaled 



norms (3.12) and (3.13) 



Theorem 3.9 (Reliability). There exists a constant C > independent of h and 
the parameters t and a such that 



\U - Uh 



■,tJi 



l\p-Ph\\ 



■,t,h 



< Ct]. 



(3.35) 



Theorem 3.10 (Efficiency). There exists a constant C > independent of h, t, 
and a such that 

<c{\\u-u,\\l,j, + \lp-pl\ll,^^ (3.36) 
+ E ( ,2.2^ ,2 11/ - MIo./. + it' + ^'hl)\\g - P,g||g,^-) }. 

Thus for the displacement Uh and the postprocessed pressure pjj we have by 
Theorems 13.91 and 13.101 a reliable and efficient indicator for an elementwise constant 
permeability parameter a and all values of the effective viscosity parameter t. In 
addition, we have the localized lower bound 



r?^ + '7l<c{||w-w,.||2,_,^^^, + |b-p;j||2_,,,,„^ (3.37) 



E ' 



in which lok C is the patch of elements surrounding an element K. and the sub- 
scripted norms above are evaluated only over the elements in lok- Thus we have a 
strong indication of the applicability of the estimator to adaptive refinement. 

4. Hybridization. A well-known method for dealing with the indefinite system 
resulting from the Darcy equation is the hybridization technique introduced in [5J 
[S]. The idea is to enforce the tangential continuity via Lagrange multipliers chosen 
suitably and relaxing the continuity requirement on the finite element space. Thus, 
we drop the normal continuity requirement in the spaces V^dm ^^^^ V/^^ and denote 
these discontinuous counterparts by V/j. In addition, we introduce the corresponding 
multiplier spaces 

Af = {A e [L\£^)r-^ I A e Pk{E), E e Eh, \\e = 0, E c dn}, (4.1) 
Af^ = {A e [L\£h)r-^ I A e Pk-i{E), E e £h, Me = 0, E c On}, (4.2) 
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in which Sh denotes the collection of all faces of the mesh. It can be easily shown, 
that the normal continuity of a discrete flux i6/j e V/j is equivalent to the requirement 



KeKh 



(4.3) 



Accordingly, the original finite element problem (3.111 can be hybridized in the fol- 
lowing form: Find {uh,Ph, ^h) G V/i x Qh x such that 



Bh{uh,Ph;v,q) + ^ {v-n,\h)aK = {f,v) 

KeKh 



(9,1), 



(4.4) 
(4.5) 



for all (v,q,fi) G V/i x Q^, x A;i. Due to (4.3), the solution {uh,Ph) of the hybridized 
system coincides with that of the original system. Thus, we need not modify the 
postprocessing procedure even if we drop the continuity requirement from the velocity 
space. 

4.1. Hybridization of the Nitsche term. However, now the matrix block 
corresponding to the bilinear form Bh{uh,Ph] v,q) is Si block diagonal system only for 
the special case t = 0, and for a non-zero effective viscosity we cannot eliminate the 
variables locally. To alleviate this problem we introduce a second hybrid variable for 
the Nitsche term of the velocity, see e.g. [14]. Recall, that the velocity- velocity term 
of the bilinear form Bh{uh,Ph', v,q) is 



aii{u, v) = (cr^u, v) + 



Eli 



MeKh 

nil 

1,K1)^^-(%|},I' 



(4.6) 



To this end, we follow [HT, and formally introduce the mean value of u-r as a new 
variable, m = ^{ui^t + M2,t-)- Thus we can write the tangential jump as 

Kl = 2(Mi,^-m) = -2(w2,^-m). (4.7) 

Now using the new hybrid variables the bilinear form ah{u, v) can be rewritten as 

afi{u,'m;v,r) — {a^u,v) ^ {(Vm, Vf) 



2a 



IK 



r)aK 



idu 



KeKh 
r)dK 



dn 



U-r - m)oK}- 



Here, the hybrid variable m belongs to a space A^/i C [L'^ (£/,,)]", the choice of which 
will be discussed subsequently. In addition, we introduce a slightly modified version 
of the norm (3.12) to encompass both the velocity and the hybrid variable: 



\{u,m)\ 



2 

a,tJi 



E 

KeKh 



a \\u 



2 

\o,K 



E I 

iKeKh 



Vm||o,k 



E 

EeSh 



^1 
hE^ 



u-r — m 



2 



(4.1 
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Since for the exact solution the jumps disappear, the biHnear form is consistent. 
Using exactly the same arguments as those presented in [TS] for (3.171, we have 

Lemma 4.1. The hybridized bilinear form ah{- ,■;•,■) is coercive in the discrete 
spaces Vfi x Mh, that is there exists a positive constant C such that 

ah{v,rn■v,rn)>C\\{v,rn)\\l^ ,^, V(t>, m) e x Al,,. (4.9) 

Note, that the stability holds for any choice of the space M.h- For complicated 
problems, this gives great flexibility. Thus, due to consistency and stability, we get 
optimal convergence rate as long as the space A^^ is rich enough. Here we choose 

Mh = {m e [L\Sh)T I Q{E)m\E e [Pk{E)]^-\, VE e £h}, (4.10) 

in which Q{E) is the coordinate transformation matrix from the global n-dimensional 
coordinate system to the local (n — l)-dimensional coordinate system on the face E. 
Let Ph : [L'^{E)]'^~^ -> Mh be the projection on the faces. We then get the 
following interpolation estimate. 

Lemma 4.2. Let u be such that u\k G [ff''+^(i4r)]" for \ < s <k. Then it holds 

\\{u~ RhU.u^ - PhU^)\\a,t,h < C{<jh + t)h'\\u\\,+i. (4.11) 



Proof. We proceed by direct computation. Scaling and the Bramble-Hilbert 
lemma [6 yield 



\{u-RhU,U-r-PhUT)\\cr,t.h 
1 



E 



< J2 '^''\\'>^~Rhu\\iii+t^ 

KdKh 

1 



\\V{u~Rhu)\\lj, 

.KeKh 



EeSh 
<c( a^h^'+^ 



„(m - J?,,m)^||o + — ||(m^ - PhU^)\\o,E 
tlE flE 



K&Kh ) 



The result is immediate after taking the square root. □ 

Combining the interpolation results with the consistency and ellipticity properties 
yields an optimal convergence rate for the velocity. 

Theorem 4.3. Assuming sufficient regularity, for the finite element solution 
(uhTirih) of the hybridized system it holds 



\\{u - Uh,u-r - mh))\\a,t,h < C{ah + t)h''\\u\\s+i 



(4.12) 



Thus the hybridization preserves the convergence rates of Theorem |3.8[ 

The residual a posteriori estimator of Section |3.4|can be modified to handle the 
hybrid variable by changing the edgewise estimator (3.331 to 



J.2 u 

2 II II 2 I E I 

^E = j^\\n,,.-mh\\o,E + ^^j;^,\ 



dn 



0,E 



a^hl + 1^ ' 



-'hlWO.E- 



(4.13) 



Following the lines of the proof in |27] it is easy to prove that also the modified esti- 
mator is both sharp and reliable. This will be demonstrated in numerical experiments 
in Section [5l 
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4.2. Implementation and local condensation. In practice, it is beneficial 
to choose the hybrid variable m slightly difl:erently, namely as the weighted average 
m — |(tti.T- + 1*2. T-). Now the hybridized bilinear form can be written as 



ahiu,m;v,r) = (a'^u^v) + V -—{m,r)gK 

E du dv 2a , . 2a , , , 

{{^.r}gK + {-^,m}oK - 7—{Ui-,r)oK - 7—{v-r,m}oK} 
on on hp, hp 

KeKh 

9 r ^ \ 2q: du du 

+ t y {{Vu,Vv)k + 7—{U-r,V^}QK - {^,V-r)dK - {^,U-r)dK}- 

,f-t hs on On 

Note, that now we get a t-independent part for the hybrid variable and the system 
remains solvable even in the limit t 0. It is clear that all of the results in Section ITTI 
hold also for the scaled hybrid variable. 

The main motivation for the hybridization procedure is to break all connections 
in the original saddlepoint system to allow for local elimination of the velocity and 
pressure variables at the element level. After hydridization the matrix system gets 
the following form where A is the matrix corresponding to the bilinear form a/i(- , • ), 
B to 6( - , • ), whilst C and D represent the connecting blocks for the hybrid variables 
for normal continuity and the Nitsche terms, respectively. AI is the mass matrix for 
the Nitsche hybrid variable. 



(A 


B^ 






B 











C 



















mJ 



(4.14) 



Since A and B are now block diagonal matrices, they can be inverted on the 
element level and we get the following symmetric and positive definite system for the 
hybrid variables. We denote the by H the following matrix that can be computed 
elementwise. 



H A-^B^{BA-^B'^)-^BA-^ - A'K (4.15) 
The matrix H is positive definite and symmetric. Eliminating the velocity and pres- 



sure from the system matrix (4.141 yields the following system matrix for the hybrid 
variables (A, m) corresponding to the normal continuity and tangential jumps, respec- 
tively. 



fCHC^ CHD^ 
\DHC^ DHD'^ + M 



(4.16) 



Evidently the resulting system is symmetric and positive definite. Note, that whilst 
the connecting block D will vanish as i — > 0, the M block does not depend on t, thus 
the whole system remains invertible even in the limit. 

4.3. Application to domain decomposition. The hybridized formulation is 
well-suited to solving large problems with the domain decomposition method. The 
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hybrid variables readily form a discretization for the skeleton of the domain decom- 
position method for any choice of non-overlapping blocks. The local problems are of 
the Dirichlet type, and the domain decomposition method can be implemented easily 
using local solvers on the subdomains. 

Furthermore, only the skeleton of the domain decomposition mesh can be hy- 
bridized using conventional i/(div)-elements for the saddle point system in the sub- 
domains. We also have great flexibility in the choice of the hybridized variables, thus 
allowing one to use a lower number of degrees of freedom on the skeleton when com- 
putational resources are limited or alternatively construct a coarser approximation 
for use as a preconditioner. Most importantly, the mass conservation property is re- 
tained over the subdomain boundaries, thus making the method a viable alternative 
for multiphase flow computations. 

5. Numerical tests. In this section, we shall numerically demonstrate the per- 
formance of the method. In all of the convergence tests we are mainly interested in 
the effect of the ratio between the Stokes and Darcy terms —t'^Au and a^u. Thus for 
the convergence tests we reduce the problem to a one-parameter family of problems 
by choosing cr = 1 on the whole of the domain fl. This allows us to visualize differ- 
ent behavior of the numerical method in the Stokes and Darcy regimes for varying 
parameter values. 

This choice of parameterization is made for clarity of the numerical analysis, since 
it allows the Brinkman problem to be easily interpreted as a singular perturbation 
of the Darcy problem, on which the flnite element scheme proposed is based. More 
exactly, the underlying Darcy flow problem for which the behaviour of the H{div)- 
conforming discretization is well-understood is kept constant as regards the ratio of the 
permeability parameter to the mesh size. The effects of adding the viscous term are 
then studied by gradually increasing value of the t parameter, which is in accordance 
with the numerical nature of the problem. 

Since the main value of the convergence tests lies in showing the performance 
of the method throughout the parameter range in an idealized case with a known 
solution, the permeability can be chosen not to vary spatially without losing generality. 
Thus, the conversion to a one-parameter family can be interpreted also as a constant 



scaling of the pressure and the loading in the original model (2.4 1 as jH] 



tV . 1 _ 1 



Au + u+ ^\7p = —f, inn, (5.1) 
div u ^ g, in Cl. (5.2) 



Accordingly, in the physically more realistic case of approximately constant viscosity 
and varying permeability, the converge results can be interpreted by regarding the 
above scaling of the pressure as a variation in the magnitude of the driving force 
exerted by a given pressure gradient. 

First, we test the optimal convergence properties of the method along with the 
performance of the a posteriori estimator in two cases with different regularity proper- 
ties. The proposed method is also compared with the Stokes- type approach using the 
MINI element. We proceed to demonstrate the effect and the importance of the post- 
processing procedure. Next, the convergence of the hybridized method is studied in 
the framework of domain decomposition methods. Our fourth test is a flow in a chan- 
nel, demonstrating the performance of Nitsche's method in assigning the boundary 
conditions and the applicability of the error indicator to adaptive mesh reflnement. 
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We end the section with a reahstic example employing permeability data from the 
SPEIO dataset. In all of the test cases the ratio h/t is the ratio l/{t\/N). in which 
N is the number of degrees of freedom. For a uniform mesh we have h/t ^ 1/ {t\/N). 
Note, that this holds only in the two-dimensional case considered in the computations. 
The data approximation term in the a posteriori estimator is neglected in the compu- 
tations. In all but the last test case in which real-life values are used, the parameters 
as well as the size of the domain are dimensionless. 

5.1. Convergence tests. For the purpose of testing the convergence rate, we 
pick a pressure p such that Vp is a harmonic function. Thus, u = Vp is a solution of 
the problem for every t > 0. In polar coordinates (r, G) the pressure is chosen as 

p{r,Q)=r^sm{pe) + C. (5.3) 

The constant C is chosen such that the pressure will have a zero mean value. Moreover, 
we have p £ H^^/^{n), and subsequently u € [H^{Q)]"', see [TS]. In the following, 
we have tested the convergence with a wide range of different parameter values, and 
the results are plotted with respect to the ratio of the viscosity parameter t to the 
mesh size h. Our aim is to demonstrate numerically, that the change in the nature 
of the problem indeed occurs ai t = h, and that the convergence rates are optimal in 
both of the limiting cases. First we choose /3 = 3.1 to test the convergence rates in 
a smooth situation. With the first-order BDM element we are expecting converge 
in the Darcy end of the parameter range and h in the Stokes limit. 

As is visible from the results in Figures [l] and [2] the behaviour of the problem 
changes numerically when t = h. Thus, even though in practical applications almost 
always i > 0, numerically the problem behaves like the Darcy problem when t < h. 
As can be seen from Figure |2] the converge rates follow quite closely those given by 
the theory. Furthermore, both the actual error and the a posteriori indicator behave 
in a similar manner. Notice, that the convergence rate exhibits a slight dip at the 
point in which the nature of the problem changes. This is a result of the dominating 
error component changing from pressure to velocity error as we pass into the Stokes 
regime. 

To show the applicability of the a posteriori indicator to mesh refinement, we con- 
sider the more irregular case f3 — 1.52. Our refinement strategy consists of refining 
those elements in which the error exceeds 50 percent of the average value. The tresh- 
old is halved until at least five percent of the elements have been refined. The edge 
estimators are shared evenly between the neighbouring elements. When compared to 
uniform refinement in Figure [4] Figure [6] shows that the convergence is considerably 
faster with adaptive refinement. Furthermore, adaptive refinement seems to allevi- 
ate the problem of convergence rate drop at the numerical turning point t — h, as 
demonstrated in Figures [3] and [5j Clearly these results indicate that the a posteriori 
indicator proposed gives reasonable local bounds and can thus be used in adaptive 
mesh refinement. 

5.2. Comparison with the MINI element. A common choice for solving the 
Stokes problem is the classical MINI element [4 . This element has been applied to 
the Brinkman problem and thoroughly analyzed both theoretically and numerically 
in |24| 1^ . We use the same test case as above with the regularity parameter set to 



(3 = 3.1. Notice, that the mesh-dependent norms (3.121 and (3.131 reduce to the ones 
presented in |24] when a continuous velocity-pressure pair is inserted. Thus we can 
use the same mesh-dependent norm for computing the error for both of the elements 
and the results are comparable with one another. 
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As can be seen from the results in Figures [7]and[8l the convergence rate for the 
MINI element is as expected of the order h throughout the parameter regime. For 
the BDMl element, on the other hand, we get convergence in the Darcy regime, 
and after a slight dip at the turnaround point convergence relative to h. Turning our 
attention to Figure [7] reveals that the behaviour of the absolute value of the relative 
error differs vastly for the two elements. Clearly, the BDM element outperforms the 
MINI element in the Darcy regime by several decades of magnitude, whereas in the 
Stokes regime the performance of the MINI element is superior. This implicates that 
it is impossible to clearly tell which element is superior for the Brinkman flow. How- 
ever, usually real-life reservoirs consist mostly of porous stone with scattered vugs 
and cracks. Thus the volume of the Stokes-type regime is often small compared to 
the Darcy regime, implying that the good performance of the BDM element in its 
natural operation conditions might offer significant performance increase for the over- 
all simulation. In problems with a greater fraction of void space, such as filters with 
large free-fiow areas separated by permeable thin layers, the Stokes-based elements 
might be a more natural choice. 

5.3. Postprocessing. In this section, we show the necessity of the postprocess- 
ing by comparing the behaviour of both the exact error and the a posteriori estimator 
for the original and the postprocessed pressure. We shall use the same test case as in 
the previous sections. First we choose f3 = 3.1 for testing the effect of the postpro- 
cessign on the convergence. On the second run we choose /? = 1.52 and use the same 
mesh refinement strategy as before to show the necessity of the post-processing for 
the usefulness of the error estimator. 

From the results of Figure |9] it is immediately evident that in the case of a small 
parameter t corresponding to a Darcy-type porous fiow the postprocessing procedure 
is of vital importance for the method to work. However, as the viscosity increases the 
weighting of the pressure error changes, and the norm is more tolerable of variations 
in the pressure. Once again, this change in behaviour appears exactly at t = h. As 
regards convergence rate, the non-postprocessed method is able to attain a close-to- 
optimal rate in the Stokes regime, cf. Figure [TO] The same pattern is evident also 
with the more irregular test case with f3 = 1.52 as shown in Figures [TT] and [T2j When 
in the Darcy regime, the indicator simply does not work in adaptive refinement since 
the pressure solution lacks the necessary extra accuracy. However, when crossing into 
the Stokes regime, convergence starts to occur, and the adaptive procedure achieves 
a rather high rate of convergence. Evidently, postprocessing is vital for the method 
in the Darcy regime. Even though the method seems to work without postprocessing 
in the Stokes regime, one cannot guarantee convergence and thus the method should 
always be used only in conjunction with the postprocessing scheme for the pressure. 
The cost of solving the postprocessed pressure is negligible compared to the total 
workload since the procedure is performed elementwise. 

Note that with postprocessing using the BDM family of elements is more eco- 
nomical than using the RT family with respect to the number of degrees of freedom, 
since we can use initially one order lower approximation for the pressure, and still get 
the same order of polynomial approximation and convergence after the postprocessing 
procedure. 

5.4. Hybridized method. Here we study the convergence of the fully hy- 



bridized method for different parameter values using the modified norm (4.8) and 



the a posteriori estimator (4.13). We employ once again the same exact solution as 
in the other convergence tests with the same values for the regularity parameter /3. 
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The errors are measured using the modified norm (4.8 1. 

First we hybridize all of the edges in 8h ■ As Figures [13] to [16] clearly show, 
the hybridized method behaves in an identical manner compared to the standard 
formulation, both in the case of a regular and an irregular exact solution. Thus it 
can be concluded that the error induced by hybridizing the jump terms inexactly is 
negligible, and the proposed convergence rates are attained. 

In the second test we divide the domain U, into 16 triangular subdomains, and hy- 
bridize the finite element spaces only on the domain boundaries and employ standard 
_ff(dz'y)-conforming BDMl elements in the subdomains. Accordingly, the estimator 
is modified only on the hybridized faces. Note, that since the hybridization for the 
tangential jumps is not exact, the system of equations for the domain decomposition 
method differs from the one for the fully hybridized method, since we are using the 



original bilinear form (3.9 1 in the subdomains. As is evident from Figures 17 to 20 
we have the correct convergence for the hybridized domain decomposition method 
too. 



Lastly, we investigate the condition number of the Schur complement matrix (4.16 1 
for /3 = 3.1 and values of t ranging from t = to t = 10^. The total number of de- 
grees of freedom is kept approximately constant, whilst the number of the subdomains 
in the domain decomposition method is varied. The subdomains are approximately 
equal in size. As is evident from Figure [21] in the Darcy regime the condition number 
is rather insensitive to the value of the parameter i, however as t is increased and 
we pass to the Stokes regime the condition number grows as t^. Furthermore, in the 
Darcy regime we observe an increasing condition number related to the size of the 
subdomain Vti as diam(r2i) ^, as one would expect. 

5.5. Pressure driven flow in a channel. A viscous flow in a narrow channel 
driven by a pressure drop gives us a test case for which the exact solution is known. 
The flow is driven by a linear pressure drop with no-slip boundary conditions for t > 0. 
For the Darcy case t = 0, only the normal component of the velocity vanishes on the 
boundary. We will test the convergence with Nitsche's method for the tangential 
boundary condition with adaptive mesh reflnement. Denote by £h,u^ the collection of 
edges residing on the part of the boundary C d^l on which we set the tangential 
velocity. To set the tangential velocity we then add the term 

Yl {i^{'^t,v^)e- {^,v^)e~ {^■,u^)e} (5.4) 
hE on on 

to the bilinear form a/i(- , • ). The loading is augmented by the term 

{T-{'iJ'D,T,V^)E - {^■,ud.t)e}- (5.5) 
p -T^ He On 

in which ud denotes the velocity boundary condition. Thus, in the limit < = a 
standard Darcy problem with Dirichlet normal boundary conditions is solved. Finally, 



the error estimator (3.33) is modifled by adding the corresponding boundary face 
terms to the estimator. 

In the unit square we take the pressure to be p = —x -\- |. Then zero boundary 
conditions for the velocity for y = and y — 1 give the exact solution u — {u, 0), with 
the x-directional velocity given by 

_ /(l + ei/*-e(i-?^)/*-e''/*)/(l + eV*), i>0 

"-\i, t = o. 
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As Figures [22] through [26] demonstrate, the adaptive process is able to catch the 
boundary layer effectively in the case where no-slip boundary conditions are desired on 
the boundary, leading to nearly identical converge rates for different parameter values 
as indicated by Figure [27] Note, that as the mesh is refined on the edges, the problem 
changes numerically to a Stokes-type problem near the boundary since the mesh size 
h drops locally below the parameter t. Different boundary conditions between porous 
media and free-flow regions are discussed from a more physical viewpoint in e.g. f l8] . 

5.6. Example with realistic material parameters. In this final section we 
consider the Society of Petroleum Engineers test case SPEIO [13] with realistic porosity 
and permeability data for a heteregeneous oil reservoir. Following [34', we make the 
common choice fl = fj,. We consider a single layer flow as a two-dimensional flow 
problem. For the outflow quantities to have meaningful units, we assume a thickness 
of 2 ft for the layer. The dimensions of the problem are 2200 x 1200 ft, the fluid is 
chosen as e.g. water with a viscosity of 1 cP. The flow is driven by a pressure on 
the left-hand side of the domain. The top and bottom boundaries have a no-flow 
boundary condition. To demonstrate the effect of the Brinkman term to the flow, we 
modify the permeability data by adding a rectangular streak of very high permeability 
with the dimensions 1100 x 20 ft in the middle of the domain. The advantage of the 
Brinkman model is the ability to model cracks or vugs by simply assigning a very 
high (or inflnite) permeability to these parts of the domain. 

In the numerical tests the performance of the a posteriori estimator for adaptive 
refinement with extremely heterogenous data with vast differences in magnitude is 
studied. We start all of the computations from a coarse triangulation of the domain 
with 328 triangles and 1352 total degrees of freedom. This results in an initial element 
diameter of approximately 150 ft. In our tests we employ layer 68 from the SPEIO 



dataset, and assume the permeability to be as in equation (2.3 1. 

In the first test case we consider the original dataset. In addition, we add a 
streak with a very high permeability of 10^ Darcy representing e.g. a gravel-filled 
crack inside the domain in two slightly different orientations, cf. Figure ]28j We also 
compute the net flow exiting the domain through the right boundary. We choose a 
simple reflnement strategy in which we always refine one percent of the elements in 
which the estimator attains the largest values. The edgewise indicators are divided 
evenly to neighbouring elements. The refinements are performed until the limit of 
100 000 degrees of freedom is reached. Figure ]29[ demonstrates that the a posteriori 
estimator is able to capture the fine details of the non-modified permeability field. In 
Figures [30[ and [3T] it is clearly visible, how the method is able to find the non-piercing 
permeability streak, even if the streak is outside the natural fiow path. In particular, 
in Figure [30[ it is clearly visible how the fiow paths to and from the high permeability 



region are found by the refinement procedure. In Figure 32 we compare the fiow rate 
for the different permeability conflgurations. As can be seen, the flow rate is stabilized 
as the mesh reflnement procedure proceeds. 

In the second test we modify the permeability field by adding a through-domain 
crack 20 feet wide into the domain by extending the streak from the first test case 
and we reduce the pressure loading to 1 Pa. The permeability is set to 10^^ Darcy in 
the crack - thus modelling in essence a void space. We compare the Brinkman and 
Darcy models. In this case we use a modified refinement strategy to speed up the 
adaptive process. In the first three refinements we refine 15 percent of the elements, 
in the next three 10 percent, followed by 5 percent in the next three and 2 percent in 
all of the subsequent refinements. We stop refining once 100 000 degrees of freedom 
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are attained. Clearly Figures 34 and 35 indicate that both for the Darcy (i.e. t — 0) 
and Brinkman models the indicator finds the piercing crack. However, as can be seen 
from the net flow results on Figure |33| the Darcy model overestimates the flow by 
many orders of magnitude due to the absence of the viscous forces. Also the velocity 
of the fluid grows unnaturally large for the Darcy model, but in the Brinkman model 
the velocity stays reasonable. The net flow rate results also clearly indicate how 
the reflnement proceeds gradually until the piercing crack is found and the flow rate 
jumps, after which it is again stabilized. 

Overall, simulating flne cracks requires extremely flne mesh around the cracks, as 
demonstrated by the mesh density plots. Thus uniform reflnement would yield systems 
too large to practically solve, in particular in 3D. We also investigated the value of the 
ratio ^Jf^^ for all of the Brinkman simulations. When this ratio is greater than unity, 
we are numerically in the Stokes regime. Clearly, in the case of high permeability 
cracks we obtain flow domains in which we actually move numerically into the Stokes 
regime in the Brinkman model when modelling fractures and cracks. 

6. Conclusions. We were able to numerically demonstrate the theoretical re- 
sults for the Darcy-based method of [27, for solving the Brinkman equation. Further- 
more a hybridization technique was introduced for the whole system, which might 
prove useful for handling very large systems with the domain decomposition method. 
The hybridized method was also shown, both theoretically and numerically, to possess 
the same convergence properties as the original problem for all values of the parameter 
t. We also demonstrated the performance of the a posteriori estimator by applying it 
successfully to adaptive mesh reflnement, and compared the Brinkman model to the 
Darcy model in the framework of an oil reservoir simulation. 
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Figure 1. Relative error in the mesh 
dependent norm. Uniform refinement for a 
smooth solution (5 = 3.1. 
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Figure 2. Converge rate for different 
values oft. Uniform refinement for a smooth 
solution f} = 3.1. 
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Figure 3. Relative error in the mesh 
dependent norm. Uniform refinement for an 
irregular solution fi = 1.52. 



Figure 4. Converge rate for different 
values oft. Uniform refinement for an irreg- 
ular solution 13 = 1.52. 
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Figure 5. Relative error in the mesh 
dependent. Adaptive refinement for an irreg- 
ular solution (5 = 1.52. 



Figure 6. Converge rate for different 
values oft. Adaptive refinement for an irreg- 
ular solution /3 = 1.52. 
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Figure 7. Relative error in the mesh 
dependent norm. Uniform refinement for a 
smooth solution f} = 3.1. 



Figure 8. Converge rate for different 
values oft. Uniform refinement for a smooth 
solution P = 3.1. 
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Figure 9. Relative error in the mesh 
dependent norm without postprocessing. Uni- 
form refinement for /3 = 3.1. 



Figure 10. Converge rate for different 
values of t without postprocessing. Uniform 
refinement for /3 = 3.1. 




Figure 11. Relative error in the 
mesh dependent norm without postprocess- 
ing. Adaptive refinement for (5 = 1.52. 



Figure 12. Converge rate for different 
values of t without postprocessing. Adaptive 
refinement for /3 = 1.52. 
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Figure 13. Relative error in the mesh Figure 14. Converge rate for different 

dependent norm for the hybridized method values of t for the hybridized method with /3 = 
with uniform refinement and /3 = 3.1. 3.1. 
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Figure 15. Relative error in the mesh 
dependent norm for the hybridized method 
with uniform refinement and /3 = 1.52. 



Figure 16. Converge rate for different 
values oft for the hybridized method with fi = 
1.52. 




r 


Estimator 




Error 


1/ 


1 




1 

1 
f 




ft 

'•~\- /I 

X /■' 
V. // 
\>-^/ 



10 

In/t ratio 



Figure 17. Relative error in the mesh Figure 18. Converge rate for different 

dependent norm for the domain decomposi- values of t for the domain decomposition with 
tion with 16 subdomains and (5 = 3.1. 16 subdomains and f} = 3.1. 
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Figure 19. Relative error in the mesh 
dependent norm for the domain decomposi- 
tion with 16 subdomains and f} = 1.52. 



Figure 20. Converge rate for 
values of t for the domain decomposition with 
16 subdomains and (3 = 1.52. 
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Figure 21. Condition number for several subdomain divisions for different values oft. 
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Figure 22. Final mesh after adaptive 
refinement t = 5 Figure 25. Final mesh after adaptive 

refinement, t = 0.05 
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Figure 23. Final mesh after adaptive 
refinement, t = 0.2 




Figure 26. Final mesh after adaptive 
refinement, t = 0.005 
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Figure 24. Final mesh after adaptive 
refinement, t = 0.1 



Figure 27. Convergence rates of the 
adaptive solution for different values oft. 
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Figure 28. Original and modified per- 
meabilities on a logarithmic scale in mD for 
layer 68 of SPEIO. 
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Figure 29. Flow in the domain with no 
streak. Original data from SPEIO dataset, 
layer 68. 
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Figure 32. Net flow rates for no streak or one internal streak. 
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Figure 33. Net flow rates for the Brinkman and Darcy models with a piercing crack with a 
permeability of 10^^' Darcy. 
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